{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import itertools\n",
    "import mdtraj as md\n",
    "import mdtraj.testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load up some example data. This is a little 28 residue peptide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t = md.load_pdb('http://www.rcsb.org/pdb/files/2EQQ.pdb')\n",
    "print(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`md.baker_hubbard` identifies hydrogen bonds baced on cutoffs\n",
    "for the Donor-H...Acceptor distance and angle. The criterion employed\n",
    "is $\\theta > 120$ and $r_\\text{H...Acceptor} < 2.5 A$ in\n",
    "at least 10% of the trajectory. The return value is a list of the \n",
    "indices of the atoms (donor, h, acceptor) that satisfy this criteria."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hbonds = md.baker_hubbard(t, periodic=False)\n",
    "label = lambda hbond : '%s -- %s' % (t.topology.atom(hbond[0]), t.topology.atom(hbond[2]))\n",
    "for hbond in hbonds:\n",
    "    print(label(hbond))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the actual distances between the donors and acceptors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "da_distances = md.compute_distances(t, hbonds[:, [0,2]], periodic=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot a histogram for a few of them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "color = itertools.cycle(['r', 'b', 'gold'])\n",
    "for i in [2, 3, 4]:\n",
    "    plt.hist(da_distances[:, i], color=next(color), label=label(hbonds[i]), alpha=0.5)\n",
    "plt.legend()\n",
    "plt.ylabel('Freq');\n",
    "plt.xlabel('Donor-acceptor distance [nm]')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
